subroutine set_xyz_grids
  use my_mpi
  use logs
  use parameter_module
  use camera_module
  use grid
  use header, only : BoxSize
  implicit none
  integer*4 :: i
  integer*4 :: NRemain
  integer*4 :: Noff

  write(logunit,*)'  -establishing smoothing grids'
  call flush(logunit)

  ! Set the number of slices on each processor
  NZLocal  = Params%nz / NumPEs
  NRemain  = NZLocal*(NumPEs-1)
  NZMaster = Params%nz - NRemain
  if(MyPE==0) NZLocal = NZMaster

  ! Set cell sizes
  CellSize(1)    = L(1)/real(Params%nx)
#ifdef CUBE
  CellSize(2)    = L(2)/real(Params%ny)
  CellSize(3)    = L(3)/real(Params%nz)
#else
  CellSize(2)    = L(3)/real(Params%nz) ! L(3) because still consider a cube
#endif

  ! Set the dimensions of the local smoothing grid
  GridExtent(1,1) = x0 - 0.5*L(1)
  GridExtent(2,1) = x0 + 0.5*L(1)
  GridExtent(1,2) = y0 - 0.5*L(2)
  GridExtent(2,2) = y0 + 0.5*L(2)
  GridExtent(1,3) = z0 - 0.5*L(3)
  GridExtent(2,3) = z0 + 0.5*L(3)

  ! z-axis planes need process-based offset, save separately
  Noff = (NZMaster + (MyPE-1)*NZLocal)
  LocalExtent(:,1) = GridExtent(:,1)
  LocalExtent(:,2) = GridExtent(:,2)
#ifdef CUBE
  LocalExtent(1,3) = GridExtent(1,3)  + real(Noff)*CellSize(3) 
  LocalExtent(2,3) = LocalExtent(1,3) + real(NZLocal)*CellSize(3)
#else
  LocalExtent(1,3) = GridExtent(1,3)  + real(Noff)*CellSize(2) 
  LocalExtent(2,3) = LocalExtent(1,3) + real(NZLocal)*CellSize(2)
#endif

  ! Store smallest value as the pixel size, as per Volker's routine
  PixelSize = minval(CellSize)

end subroutine set_xyz_grids

! ------------------------------------------------------------------------

subroutine sph2grid(itype)
  use my_mpi
  use logs
  use camera_module, only : L,x0,y0,z0
  use parameter_module
  use particles
  use grid
  use smoothing_variables
  implicit none
  integer islice, ipart, ivar, itype, idim, ierr
  integer ixcent, ixstart, ixend, iycent, iystart, iyend, ixoff, iyoff, ix, iy, iz
  real*8 :: x_cent, y_cent, rplane, dist, weight
  real*8 :: x, y, z, dx, dy, dz, dz2, dr, dr2, h, h2
  real*8 :: kernel
  character(len=16) :: cslice,ctype,cvar,cnvar

  integer*4 :: i,j,k,this_index
  integer*4 :: nx,ny,nz
  real*8    :: pos(3)
  real*8    :: hmin
  real*8    :: xx,yy,zz
  real*8    :: xxx,yyy,zzz
  real*8    :: r,r2,u
  real*8    :: kernel_sum,wk
#ifdef CUBE
  real*8    :: dCellsize(3)
#else
  real*8    :: dCellsize(2)
#endif

  dCellSize = dble(CellSize)

#ifdef CUBE
  hmin = dble(sqrt(3.0/8.0)*pixelsize)
#else
  hmin = dble(sqrt(1.0/2.0)*pixelsize)
#endif

  write(ctype,'(I1.1)')itype
  write(cnvar,'(I)')nvar_by_type(itype)
  write(logunit,*)'  -executing SPH mapping, itype = '//trim(adjustl(ctype))//' ngrids = '//trim(adjustl(cnvar))
  write(logunit,*)'   -imposing min h_sml = ',real(hmin)
  write(logunit,*)'   -imposing max h_sml = ',maxval(L)

  ! Loop over the *local* particles
  write(logunit,*)'   -smoothing local particles, N =',NumPart_ByPE(MyPE)
  call flush(logunit)
  do ipart = 1,NumPart_ByPE(MyPE)

#ifdef NO_EOS
     if(EOS(ipart) > 0.5) cycle
#endif

     ! Extract coordinates
     pos(1) = dble(Coordinates(1,ipart))
     pos(2) = dble(Coordinates(2,ipart))       
     pos(3) = dble(Coordinates(3,ipart))
     h = max(dble(SmoothingLength(ipart)),hmin)
     h = min(dble(SmoothingLength(ipart)),dble(maxval(L)))

     ! How many cells does this particle cover?
     nx = h / dCellSize(1) + 1
#ifdef CUBE
     ny = h / dCellSize(2) + 1
     nz = h / dCellSize(3) + 1
#else
     nz = h / dCellSize(2) + 1
#endif

     ! Make pos relative to the grid
     pos(1) = pos(1) - dble(LocalExtent(1,1))
     pos(2) = pos(2) - dble(LocalExtent(1,2))
     pos(3) = pos(3) - dble(LocalExtent(1,3))


     ! If out of scope, cycle to next particle
     if(pos(1)+h < 0.0d0 .or. pos(1)-h > dble(L(1)) .or. &
        pos(2)+h < 0.0d0 .or. pos(2)-h > dble(L(2)) .or. &
        pos(3)+h < 0.0d0 .or. pos(3)-h > dble(L(3))) then
        cycle
     endif

     ! x,y central pixel of region covered by the particle on the mesh */    
     x = (floor(pos(1) / dCellSize(1)) + 0.5d0) * dCellSize(1)
#ifdef CUBE    
     y = (floor(pos(2) / dCellSize(2)) + 0.5d0) * dCellSize(2)
     z = (floor(pos(3) / dCellSize(3)) + 0.5d0) * dCellSize(3)
#endif
     z = (floor(pos(3) / dCellSize(2)) + 0.5d0) * dCellSize(2)

     ! Set square of h
     h2  = h*h  

     ! Determine kernel normalizaton */
     kernel_sum  = 0.0d0
     do dx = -nx,nx,1
#ifdef CUBE    
        do dy = -ny,ny,1
#endif
           do dz = -nz,nz,1
              
              xx = x + dx * dCellSize(1) - pos(1)
#ifdef CUBE
              yy = y + dy * dCellSize(2) - pos(2)
              zz = z + dz * dCellSize(3) - pos(3)
#else
              zz = z + dz * dCellSize(2) - pos(3)
#endif

#ifdef CUBE
              r2 = xx*xx + yy*yy + zz*zz
#else
              r2 = xx*xx + zz*zz
#endif
              if(r2 < h2)then
                 r = sqrt(r2)
                 u = r / h                 
                 if(u < 0.5d0)then
                    wk = (2.546479089470d0 + 15.278874536822d0 * (u - 1.0d0) * u * u)
                 else
                    wk = 5.092958178941d0 * (1.0d0 - u) * (1.0d0 - u) * (1.0d0 - u)
                 endif                 
                 kernel_sum = kernel_sum + wk
              endif
#ifdef CUBE
           enddo
#endif
        enddo
     enddo
     
     ! If particle contributes little to kernel, just skip it
     if(kernel_sum < 1.0d-10) cycle
     
     do dx = -nx,nx,1
#ifdef CUBE
        do dy = -ny,ny,1
#endif           
           do dz = -nz,nz,1
   
              xxx = x + dx * dCellSize(1)
#ifdef CUBE
              yyy = y + dy * dCellSize(2)
              zzz = z + dz * dCellSize(3)
#else
              zzz = z + dz * dCellSize(2)
#endif

#ifdef CUBE
              if(xxx >= 0.0d0 .and. yyy >= 0.0d0 .and. zzz >= 0.0d0)then
#else
              if(xxx >= 0.0d0 .and. zzz >= 0.0d0)then
#endif                 
                 i = xxx / dCellSize(1) + 1
#ifdef CUBE
                 j = yyy / dCellSize(2) + 1
                 k = zzz / dCellSize(3) + 1
#else
                 k = zzz / dCellSize(2) + 1
#endif              

                 if(i >= 1 .and. i <= params%nx)then
#ifdef CUBE
                    if(j >= 1 .and. j <= params%ny)then
#endif                          
                       if(k >= 1 .and. k <= NZLocal)then

                          xx = x + dx * dCellSize(1) - pos(1)
#ifdef CUBE
                          yy = y + dy * dCellSize(2) - pos(2)
                          zz = z + dz * dCellSize(3) - pos(3)
#else
                          zz = z + dz * dCellSize(2) - pos(3)
#endif

#ifdef CUBE
                          r2 = xx*xx + yy*yy + zz*zz
#else
                          r2 = xx*xx + zz*zz
#endif

                          if(r2 < h2)then
                             r = sqrt(r2)
                             u = r / h
                             if(u < 0.5)then
                                wk = (2.546479089470d0 + 15.278874536822d0 * (u - 1.0d0) * u * u)
                             else
                                wk = 5.092958178941d0 * (1.0d0 - u) * (1.0d0 - u) * (1.0d0 - u)
                             endif

#ifdef CUBE
                             Value(i,j,k) = Value(i,j,k) + Mass(ipart)*real(wk/kernel_sum)
#else
                             Value(i,k) = Value(i,k) + Mass(ipart)*real(wk/kernel_sum)
#endif


                             do ivar=1,nvar_by_type(itype)
#ifdef CUBE
                                ValueQuantity(ivar,i,j,k) = ValueQuantity(ivar,i,j,k) + &
                                     Mass(ipart)*SmoothVar(ivar,ipart)*real(wk/kernel_sum)
#else
                                ValueQuantity(ivar,i,k) = ValueQuantity(ivar,i,k) + &
                                     Mass(ipart)*SmoothVar(ivar,ipart)*real(wk/kernel_sum)
#endif
                             enddo
                             
                          endif
                       endif
#ifdef CUBE
                    endif
#endif
                 endif
#ifdef CUBE
              endif   ! Only one of these two really needed...
#else
           endif      ! ...but this keeps F90 formatting tidy in emacs
#endif

#ifdef CUBE
        enddo
#endif
     enddo
  enddo
     
enddo



  ! Now loop over the imported ghost particles
if(NumGhost_ByPE(MyPE) > 0)then 
   write(logunit,*)'   -smoothing imported particles, N =',NumGhost_ByPE(MyPE)
   call flush(logunit)
   do ipart = 1,NumGhost_ByPE(MyPE)

#ifdef NO_EOS
      if(Ghost_EOS(ipart) > 0.5) cycle
#endif

      ! Extract coordinates
      pos(1) = dble(Ghost_Coordinates(1,ipart))
      pos(2) = dble(Ghost_Coordinates(2,ipart))       
      pos(3) = dble(Ghost_Coordinates(3,ipart))
      h = max(dble(Ghost_SmoothingLength(ipart)),hmin)
      h = min(dble(Ghost_SmoothingLength(ipart)),dble(maxval(L)))

      ! How many cells does this particle cover?
      nx = h / dCellSize(1) + 1
#ifdef CUBE
      ny = h / dCellSize(2) + 1
      nz = h / dCellSize(3) + 1
#else
      nz = h / dCellSize(2) + 1
#endif

      ! Make pos relative to the grid
      pos(1) = pos(1) - dble(LocalExtent(1,1))
      pos(2) = pos(2) - dble(LocalExtent(1,2))
      pos(3) = pos(3) - dble(LocalExtent(1,3))


      ! If out of scope, cycle to next particle
      if(pos(1)+h < 0.0d0 .or. pos(1)-h > dble(L(1)) .or. &
           pos(2)+h < 0.0d0 .or. pos(2)-h > dble(L(2)) .or. &
           pos(3)+h < 0.0d0 .or. pos(3)-h > dble(L(3))) then
         cycle
      endif

      ! x,y central pixel of region covered by the particle on the mesh */    
      x = (floor(pos(1) / dCellSize(1)) + 0.5d0) * dCellSize(1)
#ifdef CUBE    
      y = (floor(pos(2) / dCellSize(2)) + 0.5d0) * dCellSize(2)
      z = (floor(pos(3) / dCellSize(3)) + 0.5d0) * dCellSize(3)
#endif
      z = (floor(pos(3) / dCellSize(2)) + 0.5d0) * dCellSize(2)

      ! Set square of h
      h2  = h*h  

      ! Determine kernel normalizaton */
      kernel_sum  = 0.0d0
      do dx = -nx,nx,1
#ifdef CUBE    
         do dy = -ny,ny,1
#endif
            do dz = -nz,nz,1

               xx = x + dx * dCellSize(1) - pos(1)
#ifdef CUBE
               yy = y + dy * dCellSize(2) - pos(2)
               zz = z + dz * dCellSize(3) - pos(3)
#else
               zz = z + dz * dCellSize(2) - pos(3)
#endif

#ifdef CUBE
               r2 = xx*xx + yy*yy + zz*zz
#else
               r2 = xx*xx + zz*zz
#endif
               if(r2 < h2)then
                  r = sqrt(r2)
                  u = r / h                 
                  if(u < 0.5d0)then
                     wk = (2.546479089470d0 + 15.278874536822d0 * (u - 1.0d0) * u * u)
                  else
                     wk = 5.092958178941d0 * (1.0d0 - u) * (1.0d0 - u) * (1.0d0 - u)
                  endif
                  kernel_sum = kernel_sum + wk
               endif
#ifdef CUBE
            enddo
#endif
         enddo
      enddo

      ! If particle contributes little to kernel, just skip it
      if(kernel_sum < 1.0d-10) cycle

      do dx = -nx,nx,1
#ifdef CUBE
         do dy = -ny,ny,1
#endif           
            do dz = -nz,nz,1

               xxx = x + dx * dCellSize(1)
#ifdef CUBE
               yyy = y + dy * dCellSize(2)
               zzz = z + dz * dCellSize(3)
#else
               zzz = z + dz * dCellSize(2)
#endif

#ifdef CUBE
               if(xxx >= 0.0d0 .and. yyy >= 0.0d0 .and. zzz >= 0.0d0)then
#else
                  if(xxx >= 0.0d0 .and. zzz >= 0.0d0)then
#endif                 
                     i = xxx / dCellSize(1) + 1
#ifdef CUBE
                     j = yyy / dCellSize(2) + 1
                     k = zzz / dCellSize(3) + 1
#else
                     k = zzz / dCellSize(2) + 1
#endif              

                     if(i >= 1 .and. i <= params%nx)then
#ifdef CUBE
                        if(j >= 1 .and. j <= params%ny)then
#endif                          
                           if(k >= 1 .and. k <= NZLocal)then

                              xx = x + dx * dCellSize(1) - pos(1)
#ifdef CUBE
                              yy = y + dy * dCellSize(2) - pos(2)
                              zz = z + dz * dCellSize(3) - pos(3)
#else
                              zz = z + dz * dCellSize(2) - pos(3)
#endif

#ifdef CUBE
                              r2 = xx*xx + yy*yy + zz*zz
#else
                              r2 = xx*xx + zz*zz
#endif

                              if(r2 < h2)then
                                 r = sqrt(r2)
                                 u = r / h
                                 if(u < 0.5)then
                                    wk = (2.546479089470d0 + 15.278874536822d0 * (u - 1.0d0) * u * u)
                                 else
                                    wk = 5.092958178941d0 * (1.0d0 - u) * (1.0d0 - u) * (1.0d0 - u)
                                 endif

#ifdef CUBE
                                 Value(i,j,k) = Value(i,j,k) + Ghost_Mass(ipart)*real(wk/kernel_sum)
#else
                                 Value(i,k)   = Value(i,k)   + Ghost_Mass(ipart)*real(wk/kernel_sum)
#endif

                                 do ivar=1,nvar_by_type(itype)
#ifdef CUBE
                                    ValueQuantity(ivar,i,j,k) = ValueQuantity(ivar,i,j,k) + &
                                         Ghost_Mass(ipart)*Ghost_SmoothVar(ivar,ipart)*real(wk/kernel_sum)
#else
                                    ValueQuantity(ivar,i,k) = ValueQuantity(ivar,i,k) + &
                                         Ghost_Mass(ipart)*Ghost_SmoothVar(ivar,ipart)*real(wk/kernel_sum)
#endif
                                 enddo

                              endif
                           endif
#ifdef CUBE
                        endif
#endif
                     endif
#ifdef CUBE
                  endif   ! Only one of these two really needed...
#else
               endif      ! ...but this keeps F90 formatting tidy in emacs
#endif

#ifdef CUBE
            enddo
#endif
         enddo
      enddo

   enddo
endif

! Normalise the quantity by the weighting
write(logunit,*)'   -normalising ValueQuantity by Value'
call flush(logunit)
#ifdef CUBE
do k=1,NZLocal
   do j=1,params%ny
      do i=1,params%nx
         if(Value(i,j,k) > 1.0e-10) then
            do ivar=1,nvar_by_type (itype) 
               ValueQuantity(ivar,i,j,k) = ValueQuantity(ivar,i,j,k) / Value(i,j,k)
            enddo
         endif
      enddo
   enddo
enddo
#else
do k=1,NZLocal
   do i=1,params%nx
      if(Value(i,k) > 1.0e-10) then
         do ivar=1,nvar_by_type (itype) 
            ValueQuantity(ivar,i,k) = ValueQuantity(ivar,i,k) / Value(i,k)
         enddo
      endif
   enddo
enddo
#endif


  write(logunit,*)'  -reducing grid min and max values'
  call flush(logunit)

  ! Find the minimum and maximum values across the WHOLE grid
  ! This is useful to write out for movie making
#ifdef CUBE

  temp_min = minval(Value(:,:,:))
  temp_max = maxval(Value(:,:,:))
  temp_sum = sum(Value(:,:,:))
  write(logunit,*)''
  write(logunit,*)'   -local min(Value) = ',temp_min
  write(logunit,*)'   -local max(Value) = ',temp_max
  write(logunit,*)'   -local sum(Value) = ',temp_sum
  call mpi_allreduce(temp_min,ValueMin,1,mpi_real,mpi_min,mpi_comm_world,ierr)
  call mpi_allreduce(temp_max,ValueMax,1,mpi_real,mpi_max,mpi_comm_world,ierr)
  call mpi_allreduce(temp_sum,ValueSum,1,mpi_real,mpi_sum,mpi_comm_world,ierr)

  do ivar=1,nvar_by_type(itype) 
     write(cvar,'(I1.1)')ivar
     write(logunit,*)'   -ivar = '//trim(adjustl(cvar))//' '//trim(var(itype,ivar))
     temp_min = minval(ValueQuantity(ivar,:,:,:))
     temp_max = maxval(ValueQuantity(ivar,:,:,:))
     temp_sum = sum(ValueQuantity(ivar,:,:,:))
     write(logunit,*)'    -local min(ValueQuantity) = ',temp_min
     write(logunit,*)'    -local max(ValueQuantity) = ',temp_max
     write(logunit,*)'    -local sum(ValueQuantity) = ',temp_sum     
     call mpi_allreduce(temp_min,ValueQuantityMin(ivar),1,mpi_real,mpi_min,mpi_comm_world,ierr)
     call mpi_allreduce(temp_max,ValueQuantityMax(ivar),1,mpi_real,mpi_max,mpi_comm_world,ierr)
     call mpi_allreduce(temp_sum,ValueQuantitySum(ivar),1,mpi_real,mpi_sum,mpi_comm_world,ierr)
  enddo

#else

  temp_min = minval(Value(:,:))
  temp_max = maxval(Value(:,:))
  temp_sum = sum(Value(:,:))
  write(logunit,*)''
  write(logunit,*)'   -local min(Value) = ',temp_min
  write(logunit,*)'   -local max(Value) = ',temp_max
  write(logunit,*)'   -local sum(Value) = ',temp_sum
  call mpi_allreduce(temp_min,ValueMin,1,mpi_real,mpi_min,mpi_comm_world,ierr)
  call mpi_allreduce(temp_max,ValueMax,1,mpi_real,mpi_max,mpi_comm_world,ierr)
  call mpi_allreduce(temp_sum,ValueSum,1,mpi_real,mpi_sum,mpi_comm_world,ierr)
  do ivar=1,nvar_by_type(itype)
     write(cvar,'(I1.1)')ivar
     write(logunit,*)'   -ivar = '//trim(adjustl(cvar))//' '//trim(var(itype,ivar))
     temp_min = minval(ValueQuantity(ivar,:,:))
     temp_max = maxval(ValueQuantity(ivar,:,:))
     temp_sum = sum(ValueQuantity(ivar,:,:))
     write(logunit,*)'   -local min(ValueQuantity) = ',temp_min
     write(logunit,*)'   -local max(ValueQuantity) = ',temp_max
     write(logunit,*)'   -local sum(ValueQuantity) = ',temp_sum 
     call mpi_allreduce(temp_min,ValueQuantityMin(ivar),1,mpi_real,mpi_min,mpi_comm_world,ierr)
     call mpi_allreduce(temp_max,ValueQuantityMax(ivar),1,mpi_real,mpi_max,mpi_comm_world,ierr)
     call mpi_allreduce(temp_sum,ValueQuantitySum(ivar),1,mpi_real,mpi_sum,mpi_comm_world,ierr)
  enddo

#endif

  write(logunit,*)''
  write(logunit,*)'   -global  min(Value) = ',ValueMin
  write(logunit,*)'   -global  max(Value) = ',ValueMax
  write(logunit,*)'   -global  sum(Value) = ',ValueSum 
  do ivar=1,nvar_by_type(itype)
     write(cvar,'(I1.1)')ivar
     write(logunit,*)'   -ivar = '//trim(adjustl(cvar))//' '//trim(var(itype,ivar))
     write(logunit,*)'    -global  min(ValueQuantity) = ',ValueQuantityMin(ivar)
     write(logunit,*)'    -global  max(ValueQuantity) = ',ValueQuantityMax(ivar)
     write(logunit,*)'    -global  sum(ValueQuantity) = ',ValueQuantitySum(ivar)
  enddo
  write(logunit,*)''
  call flush(logunit)

end subroutine sph2grid

